        subroutine pseudoinverse(Sigma,U,VT,Aaux,N_eq)

        real*8 Sigma(N_eq)
        real*8 U(N_eq,N_eq),VT(N_eq,N_eq),VTT(1000,1000)
        real*8 Aaux(N_eq,N_eq)
!         real*8 APaux(1000,1000)


        ! Verificar se A = U * Sigma * V^T:

!         r = 0.
!         do i=1,N_eq
!             do j=1,N_eq
!                 Aaux(i,j) = 0.
!                 APaux(i,j) = 0.
!                 r = r + A(i,j)**2
!             enddo
!         enddo

! !         write(*,*) 'Norma de A: ',r
! 
!         do i=1,N_eq
!             do j=1,N_eq
!                 if(i .ne. j)then
!                     Sigaux(i,j) = 0.0
!                 else
!                     Sigaux(i,j) = Sigma(i)
!                 endif
!             enddo
!         enddo
! 
!         do k=1,N_eq
!             do i=1,N_eq
!                 do j=1,N_eq
!                     APaux(i,k) = APaux(i,k) + Sigaux(i,j)* VT(j,k) 
!                 enddo
!             enddo
!         enddo
! 
! 
!         do k=1,N_eq
!             do i=1,N_eq
!                 do j=1,N_eq
!                     Aaux(i,k) = Aaux(i,k) + U(i,j)* APaux(j,k) 
!                 enddo
!             enddo
!         enddo

!         r = 0.
!         do i=1,N_eq
!             do j=1,N_eq
!                 r = r + Aaux(i,j)**2
!             enddo
!         enddo
! !         write(*,*) 'Norma de A2: ',r




        ! Calculo da pseudoinversa


!         do i=1,N_eq
!             do j=1,N_eq
!                 VTT(i,j) = VT(j,i)
!             enddo
!         enddo

	do j=1,N_eq
        	do i=1,N_eq
			if (abs(Sigma(j)) .lt. 0.0001)then
				VT(j,i) = 0.0
			else
				VT(j,i) = VT(j,i)*1.0/Sigma(j)
			endif
		enddo
        enddo




	do k=1,N_eq
		do i=1,N_eq
			Aaux(i,k) = 0.0
			do j=1,N_eq
				Aaux(i,k) = Aaux(i,k) + VT(j,i)* U(k,j) ! U(k,j) jah é a transposta de U(j,k)
			enddo
		enddo
	enddo







! 
!         do i=1,N_eq
!             do j=1,N_eq
!                 Aaux(i,j) = 0.
!                 APaux(i,j) = 0.
!             enddo
!         enddo
! 
! 
!         do k=1,N_eq
!             do i=1,N_eq
!                 do j=1,N_eq
!                     if(i.ne.j)then
!                         APaux(i,k) = APaux(i,k) + Sigaux(i,j)* U(k,j)
!                     else
!                         if(Sigaux(i,i) .lt. 0.00001)then
!                             APaux(i,k) = APaux(i,k) + 0.0*U(k,j)
!                         else
!                             APaux(i,k) = APaux(i,k) + 
!      &                          1.0/Sigaux(i,j)*U(k,j)
!                         endif
!                     endif
!                 enddo
!             enddo
!         enddo
! 
! 
!         do k=1,N_eq
!             do i=1,N_eq
!                 do j=1,N_eq
!                     Aaux(i,k) = Aaux(i,k) + VT(j,i)* APaux(j,k) !esta eh a pseudoinversa
!                 enddo
!             enddo
!         enddo


        return
        end

